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Abstract 

During early 1980s, the so-called 'escape time' method, developed to 
display the Julia sets for complex dynamical systems, was exported to 
quaternions in order to draw analogous pictures in this wider numerical 
field. Despite of the fine results in the complex plane, where all topological 
configurations of Julia sets have been successfully displayed, the 'escape 
time' method fails to render properly the non-filled-in variety of quater- 
nionic Julia sets. So their digital visualisation remained an open problem 
for several years. Both the solution for extending this old method to non- 
filled-in quaternionic Julia sets and its implementation into a program are 
explained here. 

1 Introduction 

The worldwide hype for holomorphic dynamics during the turn from 1970s to 
1980s claimed the development of graphical methods to display these systems on 
the computer screen, with emphasis to complex Julia sets J |18|. On the same 
wave, thanks to the joint work by Mandelbrot and Norton, the first pictures 
of quaternionic Julia sets came also out by applying the same method as used 
in C and known under the two definitions of ' escape time" 1 or 'level sets'. Few 
works, either in graphics terms, were published on this subject since early 1980s: 
namely, the book J3| an( A the articles by Norton Holbrook [I51EI], by 

Hart, Kauffman and Sandin [7||Sj. For analogous studies in higher dimensional 
dynamics, see Quaternions are multi-dimensional values consisting of 4 
vectors in the form 

h = r + mi + nj + pk, r,m,tj,pelR, i 2 = j 2 = k 2 = — 1 (1-1) 

and generate a non-abelian group. The quaternion space is named H for conven- 
tion, in honor of William R. Hamilton (1805-1865), inventor of such numerical 
values for solving some issues arising from rotation in 3-D Euclidean spaces. 
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Basin to oo 




Figure 1.1: Filled-in Julia sets J. Both topological and dynamical were explained 
after the first explorations of iterates in the complex numbers, by means of the 
quadratic map f(z) : z 2 + c. Orbits outside J converge to oo and they are best 
viewed in the Riemann sphere model of picture [2, 21 

The ' escape time'' method enjoys one major limitation: as one tries to display, 
for example, the non-filled-in Julia sets of basins B of attraction to a finite limit 
fixed point 8 and extending to oo (often retrieved by the iterates of rational 
maps). This article explains how to modify the 'standard' 1 method so to fill 
this leak and extend it to the 'non-fillcd-in' case. After resuming the technical 
details of the standard method in section |2 we discuss the novelties in section 
|3J finally, the implementation into the program QHD is explained in section FOI 

2 The background 

2.1 Julia sets topology 

Up to complex numbers and given a rational map 




where P, Q are co-prime polynomials, i.e. with no common roots, it was observed 
that Julia sets J (in any number field) may: 

1. be (pointwiscly) totally disconnected; 

2. have n — 1 dimensions; 

3. have n dimensions. 

QUESTION 2.1.1. Could this classification hold for any n- dimensional number 



In K., where n = 1, points 1) and 2) are the same, so that real Julia sets 1Z 
are pointwisely totally disconnected (in fact, for point 2), we have: n — 1 = 0) 
or the whole real line (for point 3), TZ is n = 1 dimensional). 

1 We liked to use this naming convention in order to differ the ' escape time' method from 
our new version. 



field? 
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In C, where n — 2, J could be either a dust of points (O-dimensional), also 
a line (n — 1 = 1 dimensional) or an area (n = 2 dimensional) . 

As far as the author knows, it is not known whether point 1) holds for 
quaternionic J too, or they are just disconnected sets with (hyper-)volumes. We 
will implicitly assume here the first two cases only, discarding a priori the third 
one, which leads to the trivial case where J spans over the whole 4-dimensional 
space. 

2.2 The open problem 

For complex numbers, any Julia set J can be successfully displayed with the 
already available methods; new and customised methods are anyway required 
to explore local vector fields 2 about the fixed points S. But quaternionic J did 
not share the same luck. The method employed in H was radically borrowed 
from C. However, as the related literature of those times indicates - together 
with the book 0, the road-map was to display the iterations of polynomials and 
to develop the details enhancement: in fact, those graphical examples always 
concerned of disconnected or of the filled-in varieties J only, mostly coming 
from the family of general quadratic polynomials 

ph 2 + q, h,p,qeM. (2.2) 
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Figure 2.2: Non-filled-in example. A case arising for 
the map IQl : this view is taken on the Riemann sphere, 
in the neighborhood of the point at infinity. 

Julia sets are mostly a mathematical question. This flaw in graphical terms 
was just the side effect from the interests of those pioneering researchers, who 
mainly wanted to start from the easier properties of polynomials 3 . Quaternions 
are non-commutative and the study of rational maps in H may get to high level 
of complication, in mathematical terms: so the investigation on these maps was 
not pursued. So it came that the graphics development was not considered 

2 Like for the so-called 'hedgehogs', invariant sets about irrationally indifferent fixed points 
not enjoying the diophantine condition or rational approximation; the digital display was also 
an open problem here. A related graphical method arranging a good solution is going to be 
explained in [2j. 

3 In the course of these early mathematical explorations, they proved that quaternionic J 
display in 3 dimensions as rotations of the related J in C around the X-axis. 
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along that direction: there was no need to quest on the central question of this 
article throughout the older works, because all topological cases of Julia sets 
for polynomials were already fulfilled by the standard method. (The complete 
related literature was listed in the bibliography.) Thus we can reasonably state 
that this problem was still open. A filled-in Julia set J includes all seed points 
whose forward iterated orbits do not tend 4 to oo. Topologically speaking, the 
boundary of a filled-in J is homeomorphic to a closed and bounded curve (even 
with fixed points). But exceptions are thrown when J also extends to oo: 
for example, in the classic example of the cubic rational map, resulting from 
Newton's method applied to the complex function f(z) : z 3 — 1, 

m : («) 

This variety of Julia sets J will be said 'non- filled-in': this case is out of the 
range of possibilities for the 'escape time' method when applied to quaternions. 

Being J displayed as a three dimensional model, the latest improvements, 
dating back to early 1990s (see references list), mainly concerned of the applica- 
tion of efficient illumination techniques or of improving the boundary refinement. 
Successful results in the former direction were achieved by means of Phong il- 
lumination model |Sj ; whereas the reader may consult J7J |H] about the second 
issue. Today software performs higher quality pictures of J than ever before, 
but these programs 5 still rely on the standard method: thus only the filled-in 
J can be properly displayed. After deeply scanning the Web and consultations 
with experts in this topic, programs displaying the non-filled-in family have been 
not obviously coded yet. Author's solution wants to overcome this limitation; 
the related code has been already implemented and tested to work into one PC 
program 'QHD', freely downloadable from author's site 6 and applied to render 
some among the further pictures. 

2.3 The standard method 

From (|1 . 1|) . one understands that Julia sets J exist in four dimensions. So 
one is not able to figure them in mind or to display somewhere. One is just 
able to embed them into two (M 2 ) or three (]R 3 ) real Euclidean spaces onto the 
computer screen. In this sense, different two-staged methods are available but 
in general they all root to a first stage where values/points are computed and 
to a second one, when the screen points are colored. In C colors palettes may 
be arbitrarily set (even to render purely aesthetic effects), but in H rigor shall 
be followed because J shall give their solid model look. 
In practice, this translates into this three steps approach: 

4 This is not the same as stating that these orbits should converge to a constant limit fixed 
point <5, due to the complicate dynamics — where one cannot properly speak of convergence - 
in the neighborhood of indifferent fixed points, both for the rational and irrational case. 

5 In fig. ??, the reader can watch the earliest results from the application of our method 
by one of such high quality programs. 

6 http:/ /malilla. supereva.it 
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Figure 2.3: A region of the 3-D Euclidean space is scanned, so 
that each point is then associated to one quaternion (by setting 
the X,Y,Z values to r, m, n respectively and p = 0). 



1. to scan a given region, provided a discretization of it in order to have a 
n— dimensional grid of values p; 

2. to iterate each p for a number of times fc; 

3. to test the resulting iterated value Pk in order to color them differently or 
choose which ones have to be plot. 




Figure 2.4: Examples of filled-in Julia sets J for the quadratic h + c. 



The 'standard' method was applied by Mandelbrot and Norton and by Hol- 
brook in their first visualizations and displays the 3-dimensional models of J . 
We will follow the same line here, focusing on Julia sets J possibly showing up 
as totally disconnected or 3-D surfaces. The care of illumination is fundamen- 
tal for any 3-D model, in order to prevent the flat looking. Phong approach 
guarantees high performance, for example, with regard to details, enhanced as 
the shining effect is rendered very accurately and very close to real light. In 
general the standard method, applied to J7, is son of what usually denoted 'es- 
cape method? in the early literature |18j devoted to the digital visualization of 
J . We pointed out that malfunctions come up for non-filled-in J from iterates 
of rational maps, for example from the application of the Newton-Raphson's 
method. In this case, one might find out those basins of attraction to a finite 
limit point and extending to oo. See p. 0] including a picture from the cubic 
quaternionic function f(h) : h 3 — 1, which turns into: 



f(h) _2h 3 + l 



(2.4) 



The Julia set curve J includes infinitely many fixed points S and also extends 
to infinity. As the standard method applies to (|2.4(l . one sees fig. Ij3.6|l . One 
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Table 2.1: Experiments with the standard method. All top pictures 
have been rendered by 'escape time' method. Columns relate to bailout ra- 
dius, whereas rows to iterates 5, 12, 24. The bottom picture was computed 
by the new 'cut-off rate' approach for watching the correct rendering. 

also notices how roughly it works in this situation: the shape of J is not so 
clear and just a part of it is displayed, but imprecisely. The grey region is the 
evidence of the failure of the standard method, which cannot adequately check 
the sequences of iterated values here. 

2.4 Tuning the parameters 

Before entering the heart chapter, one should be assured that a modification 
really urges. So we tried to tune the two parameters involved in this kind of 
experiments: the bailout disc radius and the number of iterations. Finally, we 
reported the resulting pictures in table l^TTl This experiment is relevant because, 
for sake of completeness, one shall either show the reasons why something is 
required and why it is not: these pictures want to be the empirical evidence of 
what will be further discussed in more details at section EQ1 The table ETTI at 
p. HO includes figures of the same Julia set J , the large figure at the bottom 
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is correctly drawn by the 'cut-off rate' approach and shall be compared with 
the small ones, rendered by the 'escape time': the number of iterates and the 
bailout circle radius (centered at the origin in the classical method) were tuned 
to distinct values. In some cases, one finds out a solid piece (the radius at 1.0 
was even too large and gave a great approximated shape); for larger radii, the 
grey region indicates that the location of J cannot be distinguished from the 
basin locus, because the detection rule failed 7 . Finally, the black pictures denote 
that the chosen conditions were rather insufficient. 

3 Opening quaternionic Julia sets 
to rational maps 

3.1 Disclaiming the Analysis 

After implementing the currently available methods to display quaternionic Ju- 
lia sets J for polynomials, the author wanted to achieve the same even for 
iterated rational maps. In mathematical terms, the study of the algebraic prop- 
erties related to them may be harder than in R or in C, owing to the non- 
commutativity property enjoyed by quaternion numbers. The next sections 
have been deliberately compiled to solve a practical problem exclusively, not 
regarding the mathematics behind the nature of the functions involved. In fact, 
the question on quaternionic Julia sets connectivity still deserves a theoretical 
systematisation, because the arising topological configurations are not quite the 
same as in C, for example disconnected Julia sets in H do show connected sub- 
regions, whereas they are totally disconnected (Cantor dust) in C. This also 
invalidates the same definition of Mandelbrot set as it holds in C. But this goal 
is out of scope for this article and it goes beyond author's knowledge. 
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Figure 3.5: Julia set for ITU in C. Fi gure 3.6: Malfunctions. The quater- 

nionic Julia set generated by 12.41 and pro- 
cessed by the standard method. 



7 See section [3. 31 for technical details. 
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3.2 Starting considerations 

The standard method scans the given area (in C) or volume (in R 3 , the reduced 
vector space from HI), associating a numerical value per each element of this 
region; the value is iterated for a number of times under a given function. If the 
resulting value is still inside a sphere of given radius, then the orbit is said to 
be 'trapped' and then it belongs to a basin of attraction bounded by a Julia set 
J . But this only works well when J is a closed curve which the test disc is able 
to catch completely up; on the contrary, the disc does not fit any more for the 
whole J extending to oo. Historically speaking, the standard method inherits 
the obvious limitation from the approach devoted to display Julia sets J on a 
computer and whose location and topology inside C were a priori known, either 
if filled-in or disconnected. In fact, the iterated maps, whose J were going to be 
shown onto a computer screen for the first times ever during early 1980s, were 
quadratic polynomials and namely in the form 



here all cases of J , their topology and the attracting fixed points are known 
to approximately lie inside a bailout disc with radius 2 at most. The standard 
method was invented with these 'bounded' dynamics in mind and it fits the 
range of cases offered by (|3.5|) only; thus, in this sense, it is a customised method 
which unfits the cases not covered by its possibilities. 

The standard method can be reset by changing the radius of the trapping 
disc in C but limitations come up as one deals with quaternions: such method 
strictly depends on the same functionality of the bounded ball (as the concept 
of trapping disc turns into for the embedding of H into R 3 ): i.e. when orbits are 
tested to stay bounded or to escape to oo; in fact the same method implicitly 
assumes that there is only one finite and attracting fixed point and that its 
location, together with the related basin, is a priori known to lie inside the 
test ball itself. So, for more general purposes - i.e. regardless of fixed points 
locations and quantity - the standard method does not work well. 

This will be the starting point upon which the enhancement will be based. 



z 2 + c, 



z, c G C; 



(3.5) 



No boundary 






Figure 3.7: First step. Quatcrnionic Ju- 
lia set boundary detection. 



Figure 3.8: Second step. The recur- 
rent process for approximating the Julia 
set point location. 
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3.3 Crash course reasons 



Therefore why does the standard approach fail here? 

First we say why it does not in C: actually because all points of a given 
region are iterated and tested, regardless of their location and even visibility is 
not a relevant condition; in fact one draws both the basin and the Julia set in 
C. But in 3-D, one needs to differ them because the basins shall be not drawn 
or one would just watch a cubic block instead ! The points of quaternionic J do 
need to be scanned and sorted if they are visible to be plot or not. This is the 
main direction which methods for quaternionic Julia sets must follow straight. 
The exclusion of the basins during the plot is a delicate feature. 

'Filled-in' or disconnected Julia sets configurations keep an element making 
the difference: the basin Boo to infinity. Its role is clearer in the 'filled-in' case 
(fig. Il.lf /B): one plays a boolean situation where points of B^ exit the test 
sphere, others do not. The rule is that J is met when, given two close seeds 
to J, the related orbits are checked by the 'escape time' and one finds that an 
orbit runs to infinity, the other to the finite attracting point (fig. 13. 7[) . So the 
spatial interval, wherein a point of J should be approximately located, shrinks 
over and over again up to lower bound, while the fates of the two orbits keep 
different (fig. EHJ. 




(A) (B) (C) 

Figure 3.9: The above figures resume three topological configurations for com- 
plex Julia sets: (A) totally disconnected and with the only basin to infinity, 
(B) connected and splitting C into two basins only, a basin to infinity and one 
to the finite fixed point. Then (C) a Julia set and all basins (two or more) of 
attraction to finite points. 

The standard method crashes because, for such maps as Q2.4JI . infinity is no 
longer attracting and J splits H into basins of attraction to finite 6: as J is 
even crossed, it is not detected because neighboring orbits never escape the test 
sphere, and the above rule blows away (see fig. 13. 6[) . 
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3.4 The 'cut-off rate' method 



The modification mainly relies upon the same region scanning process like in 
the standard method, but it adds the previous considerations, together with the 
inspiration from the application of normal families to holomorphic dynamics. 
Let / be an holomorphic function, then 

f(z)=z,f 1 = f(z),...f n (z)=f n - 1 (f(z)) 

are said to be the iterates of rank 0, 1, . . . , n. From iteration theory, one knows 
that, for a same attracting fixed point jk, the total basin B lk of attraction is 
the union set of all the points whose orbits converge to 7^: 

lim /™0) = 7, Vz G B lk . 




Figure 3.10: Iterates whose distance is smaller than a given value 
are not considered for the final plot. 

Let T>i , T>2 , . . . T> n , be Jordan curves bounding simply connected domains 
£>i, D2, ■ ■ ■ , D n around jk respectively (see fig. 13.100 . so that -D„-i D D n ; thus, 
in metric terms, the distance d n , of any pair of points p n -i S D n _i,p n G D n 
(points of an orbit), shrinks to as the orbit itself gets closer to 7^: 

d n = \p n —Pn-i\ < e, e>0,lime = 0. (3.6) 

So, generally speaking, given a seed point w € B lk and close to J ', the 
distance d n — > as the iterative index n — > 00. Conversely - very important to 
understand how the modification works -, one may also state that if d n is set 
to a constant s, then one can re-elaborate 13.611 into 

1/nH ~ f n +l(w)\ > S 

as the seed w gets closer and closer to J . This offers the chance to part all 
possible orbits generated throughout one basin £> 7fc into two groups: 

1. the orbits of the seed points p which converge, to the fixed point 7^., under 
a given distance d n and after a smaller number of iterations; 
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2. other orbits of seeds p which need a larger number of iterates to converge 
under d n . 

The points being closer to J belong to the second group. Hence the goal of 
revisiting the standard method is the 'isolation' of the second group of points 
and then plot them exclusively, with regard to other basins distributions where 
even the test ball of the standard method, with fixed location, does not fit any 
longer to return clear and good graphical results. If the test ball condition is: 

\Pn ~Pn-l\ < r - (3-7) 

This isolation can be simply performed by assuming the reverse inequality 
of (EH) , that is: 

\p n ~Pn-i\ > r. (3.8) 
The operator '>' achieves a reverse task: rather than trapping, external 
orbits to the ball are assumed or in other terms, the successive iterated domains 
D n exclusively, whose distances d n > r. That is, the domains of points which 
are farther, in terms of iterates, from 7^: conversely speaking, the points which 
are closer and closer to Julia sets. Moreover one drops off the test condition of 
the escape method, about the upper bound of the iterative index, which prevent 
infinite loops; we preferred to adopt the first member of inequality (|3.8|1 . because 
it showed to work more finely: the major enhancement is to offer, by means of 
just one condition of convergence, the work on n test balls, as many as the 
number of attracting fixed point 7^ of the iterated map. Therefore the reader 
has no longer to deal with the issue of knowing a priori the number and the 
location of the attracting fixed points S or how long one orbit takes to converge 
in some sufficiently close neighborhood of the fixed point. One sees such test 
balls, which do not depend on the location anymore, rising up from i|3.8[l and 
they work at the same time without their centers to be preset ad hoc in the 
algorithm code. This condition grants any arbitrary input map. 

The black regions in fig. (??) indicate that those points have not been 
deliberately plot, as result from the isolation task achieved by the 'cut-off rate' 
method. 

4 Related software 

In this section we will give a short overview on some related software which is 
currently implementing the 'cut-off rate' method. 

4.1 QHD : Quaternionic Holomorphic Dynamics 

First we show one author's program, coded to make the early experiments on 
the 'cut-off rate' method: QHD is a branch of another application, coded years 
before and named 'Inwards to Chaos', which is devoted to iterations in real, 
complex and quaternion fields. In order to focus on the visualisation techniques, 
QHD came up and became the natural environment to display the graphical 
results related to the 'cut-off rate' method. 
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Figure 3.11: Cut-off rate method. The quaternionic Julia set for 

2h z - 2h 2 - 3 
3h 2 + Ah + 1 ' 

The set was sectioned to notice either the rotational symmetry with the 
Julia set in C and the existence of interweaving basins which would make 
the classic method crash. 

4.2 The visual interface 

This application opens with a user-friendly interface to choose the set to draw 
('Julia' or 'Mandelbrot'), the drawing algorithm ('Standard' or 'Extended'), 
a default set of formulas and the illumination method ('Simple Lambertian', 
'Lambertian' and 'Phong'). 

When properties are displayed, the same window expands and shows ad- 
ditional parameters. First the details level yielding as finer pictures as it is 
incremented (as well as computation times). Then the number of iterates, the 
degree of predefined maps, the region coordinates and the value of parameter q, 
like in formula l2.2[l . 

Input features also list the command to load/save the parameters and the 
metrics related to Julia set configuration, such as the bailout value playing as 
the radius of the disc entrapping the orbits in the standard method or as the 
value r in the inequality (|3.8|) of the extended approach. 

The main interface offers a very little preview shot by a quick plot, but one 
can display magnifications into a larger window; a set of four arrows allows 
finally to change the input region coordinates. 

Additional features include a player to stop the drawing process manually, a 
window to rotate the given object along the three X-Y-Z axes, together with one 
panel to manage all parameters related to the illumination models (location of 
light points, intensity, ambient diffusion). Then, there is the possibility of setting 
any triplet of the four quaternionic components into the X-Y-Z coordinates 
system in order to watch the possible shapes assumed by quaternionic Julia sets 
into a reduced 3-D vector space (see examples in table at p. 1131) . The 
output allows to save the current picture into one graphics format file or to be 
copied into the computer clipboard for exchanging it with other applications. 
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Figure 4.12: On the left, the main interface of QHD. 

4.3 Looking around 

My goal was to test if that method really worked in other applications out- 
side mine and verify if my results were achieved in particular situations or due 
to possible mistakes in the code. In the recent period (February 2006), af- 
ter testing the method on QHD, I wanted to let the 'cut-off' method circulate 
among groups interested in computer graphics devoted to quaternionic Julia 
sets, which are usually frequented by 'mathematized' programmers, that is, 
non-mathematicians, strictly meaning, but skillful people liking to deal with 
Mathematics. One is Terry W. Gintz who welcomed the new method and im- 
plemented it into a program, retrieving high-quality renderings of such fractals; 
for example see the pictures in the table ?? at p. ??. 

5 Conclusions 

The program 'QHD' can be freely downloaded from the web-site address at 
http://www.malilla.supereva.it; Opinions and suggestions are welcomed 
to enhance it. The C++ code can be asked via author's e-mail. 

Alessandro Rosa 
Brindisi, Italy 
zandor_zz@yahoo . it 
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